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We have studied the rationally frustrated Josephson-junction array model in the square lattice 
through Monte Carlo simulations of (2 + 1)D XY-model. For frustration / = 1/4, the model at 
zero temperature shows a continuous superfiuid-insulator transition. From the measurement of the 
correlation function and the superfluid stiffness, we obtain the dynamical critical exponent z = 1.0 
(N . and the correlation length critical exponent v = 0.4 ±0.05. While the dynamical critical exponent is 

$H ' the same as that for cases / = 0, 1/2, and 1/3, the correlation length critical exponent is surprisingly 

Oh, quite different. When / = 1/5, we have the nature of a first-order transition. 

The Josephson-junction array model has attracted a great deal of attention as a prototypical model to describe 
the zero-temperature quantum phase transition jl[] in two dimensions. The competition of the charging energy and 
the Josephson coupling energy yields the transition between superfluid and an insulator. Recent development of the 
^ \ nanotechnology makes it possible to tune the transition by changing the capacitance and the Josephson coupling of 
the arrays in experiments [Q. The transition is also tuned by applying a magnetic field perpendicular to the array 
The magnetic field generates frustration in the phase coupling between the junctions. The feature of this model 
describing the quantum mechanical phase transition has been widely adopted to explain the superconductor-insulator 
transition in thin films Q|. 

The zero-temperature phase transition of the unfrustrated (/ = 0) Josephson-junction array model belongs j^, [7]] 
to the 3D XY~-model universality class with the dynamical critical exponent z — 1 and the correlation length critical 
exponent v = 0.67||. The transitions of the model with frustration / = 1/2 and / = 1/3 have the critical exponents 
numerically very close to those for / = 00. It is useful to treat the model at / = 1/q in terms of q coupled unfrustrated 
XY-Totois (i.e., order parameters) for small g[|l0[ [llj]. In terms of vortices, at / = 1/2, for example, vortices form 
two different stable structures and the transition is mostly derived by excitations with respect to them. However, the 
' couplings between rotors is believed to be always relevant. As a result, one may expect different critical exponents 
. for different value of q. 

\Q ' Furthermore, when frustration is very small so that the number of order parameters is large, the first-order nature 
of the transition will emerge. This situation is analogous to the Potts model which shows the crossover from the 
continuous transition to the first-order transition as the number of the flavors increases j^. This conjecture is 
supported by the numerical simulation of the 3D XY-mode\ with frustration / = 1/6 in a triangular lattice |l3], 
showing a first-order transition. It is also supported experimentally by the first-order melting of dilute vortex-lines in 
high-T c superconductors jl4| , whose underlying physics can be captures by the 3D JTY"-model with small frustration. 
It is, therefore, interesting to investigate the nature of the transition and to measure the associated critical exponents, 
C$ ' if it is a continuous phase transition, in the weakly frustrated Josephson-junction array model. 

In this work, we study the zero-temperature phase transition of the Josephson-junction array model in a square 
lattice when the frustration is / = 1/4 and / = 1/5. The model can be mapped into (2 + 1)D classical XY~-model with 
frustration on the xy-plane. Through Monte Carlo simulations, we find that when / = 1/4 the superfiuid-insulator 
transition of the model is continuous with the dynamical critical exponent z — 1.0 and the correlation length critical 
exponent v = 0.4 ± 0.05. While the dynamical exponent is numerically the same value as found for cases / = 0, 1/2, 
and 1/3, the correlation length critical exponent is surprisingly different. When / = 1/5, the probability distribution 
near the transition strongly suggests that the nature of the transition is of first-order. 
The Josephson junction array model is represented by the Hamiltonian 
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&\ H = Ecn 2 , -Y, E J cos (^ - °i ~ A n) (!) 

1 (ij) 

where E c and Ej denote the charging energy and the Josephson coupling energy respectively, 9i the phase of the 
superconducting order parameter at site i, rii = (l/i)(d/d9i) the number of particles at site i, the phase shift due 
to the magnetic field, and (ij) pairs of nearest neighboring sites. Frustration / is embedded in Aij whose plaquette 
sum is constant over the whole lattice, Aij — 2nf. It is convenient to use the notation i = (i x , i y ) to represent the 
position of the i-th site. Using the Landau gauge, one may have 

,4 r = /± 27r /^ if 3 = (ix,iy ± 1); 
* 3 1 otherwise. 
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FIG. 1: The probability distribution as a function of the bond energy of the classical action (a) at / = 1/4 and (b) at / = 1/5. 
The double maxima feature, the signal of the first-order transition, appears clearly in (b), but it is absent in (a). 



When E c 3> Ej, the system has the insulating ground state, while in the opposite limit it has the superfluid ground 
state. One may treat this transition more conveniently using the equivalent classical effective action]7| 

S[9] = -Kjyoafa - 9- i+5 J + cos(^ - 6~ i+Sy - 2nfi x ) + \ cos(0; - 6 1+&t )\ , (2) 

i 

where i = (i X) i y , i T ) denotes the discrete lattice sites of the (2+l)D lattice (i + S x = (i x + 1, i y , i T ), etc. ), and K the 
tuning parameter which effectively changes the ratio of the Josephson coupling energy to the charging energy. We 
multiply the coupling in the r-direction by a constant factor i to make the superfluid stiffness and the compressibility, 
which are discussed below, comparable. Here the partition function is given by 

Z = Tr e e- S[e] . (3) 

We simulate the model given by the classical action in the systems whose sizes are represented by L x L x L T . 
The first issue we are interested is the nature of the transition. In terms of vortices, the transition is described as 
the vortex-line melting. We measure the probability distribution near the expected critical point with respect to the 
bond energy, e, of the classical action, defined by 

' ^fcosfft - 6^ x ) + cos(fl ; - 6~^ x - 2wfi x ) + \ cos(^ - e~ i+Sr )] , (4) 
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Figure |l| shows the probability distribution when / = 1/4 and / = 1/5. Double maxima, separated by the distance 
of the latent heat and corresponding to the two coexisting phases in a first-order transitionjl2|, do not appear when 
/ = 1/4. This strongly suggests that we have the continuous nature of the transition in this case. When / = 1/5, 
however, the double maxima in the probability distribution do appear, indicating that the nature of the transition is 
of first-order. 

For the continuous phase transition at / = 1/4, we measure the critical exponents. Since the transition is a quantum 
phase transition, the correlation length diverges along the (imaginary) temporal direction as well as along the spatial 
direction. This requires prior knowledge of the dynamical critical exponent to determine the aspect ratio, L T /L Z , 
of the finite size systems in order to use the finite-size scaling analyses. The dynamical critical exponent can be 
determined from the measurement of the correlation function. In finite size systems, the correlation function in a 
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FIG. 2: Log-log representation of the correlation function along a spatial direction, C x (r), and the imaginary temporal 
direction, C r (r), measured in the lattice whose size is 32 x 32 x 32. The solid lines represent the fitting curves. Because of the 
frustration the spatial spin-spin correlation oscillates as a function of the distance between spins. When f — 1/4 the periodicity 
is 4, as shown in the figure. This means that only those spins separated by multiples of 4 lattice constants are correlated each 
other in phase. The fitting curve of C x (r) represents the correlation of those spins. 



spatial direction, say x, and in the imaginary temporal direction will follow the form 

C x (r)~r-"- + (L-r)- 1 '- , (5) 
C T (r) ~ r~ Vr + (L T - r)- y - , (6) 

respectively with y x = (d — 2 + z + rj) and y T = (d — 2 + z + rj)/z. Here d = 2 is the spatial dimensionality. Figure |2| 
shows the correlation function measured at K — 1.296 in the lattice of 32 x 32 x 32. Here we obtain that y x w 1.054 
and y T w 1.049, implying that z = 1.0 and rj « 0.05. 

In finite size systems the superfluid stiffness, the stiffness in a spatial direction, say x, scales with respect to the 
diverging correlation length. In the Landau gauge used above, the superfluid stiffness is given by 

i i 

where (• • •} means the average over the action of Eq. (2). The finite-size scaling ansatz of the superfluid stiffness can 
be expressed asf7| 



l —p{L^5,L T /L z ) , (8) 



P Ld+z-2 

where p is a scaling function, and 8 = K — K c is the tuning parameter. The compressibility, which is the stiffness 
in the imaginary temporal direction, is another quantity showing the scaling. The scaling behavior of compressibility 
will follow 

K = j^*(LV»5,L T /L*) , (9) 

where again k is a scaling function. 

In order to obtain the scaling behavior depending on the size of the systems, we assume z = 1 and take lattices 
of four different sizes: i.e., 20 x 20 x 20,24 x 24 x 24,28 x 28 x 28, and 32 x 32 x 32. Fig. | shows the scaling 



4 




FIG. 3: The finite-size scaling behavior of the superfluid stiffness and the compressibility. Curves for different sizes cross 
at a point, supporting that the dynamical critical exponent is z = 1.0. Inset: The curves collapse into a single curve with 
v = 0.4 ± 0.05 with respect to the scaling variable L 1 ^ (K — K c ). 



behavior of the superfluid stiffness and the compressibility. The curves of L z p and L k for many different L cross 
at K c = 1.2958 ± 0.001, implying that the chosen value of the dynamical exponent is correct. The correlation length 
exponent v is extracted by plotting the curves with respect to L X I V (K — K c ). We find that with v = 0.4 ± 0.05 the 
curves of L z p collapse onto a single scaling function quite well as shown in the inset of Fig. ||. The same critical 
exponents are obtained from the scaling behavior of the comressibility. It is interesting to note that when / = 1/4 
the value of v deviates from the one for / = 0, / = 1/2, and / = 1/3. 

In summary, we study the superfluid-insulator phase transition of the weakly frustrated Josephson junction array 
model in two dimensions. Using Monte Carlo simulations, we investigate the nature of the transition and measure the 
critical exponents from the scaling behavior near the transition. We find that when / = 1/5 the nature of transition 
is of first-order. When / = 1/4, we have the continuous transition with the dynamical critical exponent z = 1.0 and 
the correlation length critical exponent v = 0.4 ± 0.05. The value of v is different from the one for the transitions at 
/ = 0,1/2 and 1/3. 
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